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ABSTRACT 

Acoustic  sensors  offer  the  potential  to  detect,  track,  and  classify  enemy  ground  tar¬ 
gets  such  as  tanks,  trucks,  and  other  military  vehicles.  Traditional  multitarget  tracking 
techniques  partition  the  track  estimation  operation  into  two  isolated  processes:  direction- 
of-arrival  (DOA)  estimation  from  array  snapshot  data,  followed  by  track  estimation  from 
the  DOA  estimates.  In  this  paper,  a  multiple  target  track  estimation  method  that  operates 
directly  from  broadband  array  data  is  presented.  The  maximum  a-posteriori  (MAP)  esti¬ 
mator  for  contact  states  is  derived  for  temporally  uncorrelated  broadband  target  signals  and 
uncorrelated  target  tracks.  This  technique  extends  the  MAP  multitarget  tracking  approach 
developed  for  narrowband  signals  in  [1],  The  track  estimator  is  an  iterative  algorithm  em¬ 
ploying  a  nonlinear  programming  penalty  method  in  conjunction  with  an  alternating  max¬ 
imization  algorithm  for  obtaining  penalized  maximum  likelihood  (PML)  DOA  estimates. 
The  penalty  function  couples  the  DOA  estimates  from  the  PML  algorithm  to  the  tracker 
as  synthetic  measurements,  eliminating  the  data  association  step  of  traditional  multi  target 
tracking  approaches.  It  also  creates  a  feedback  mechanism  to  enhance  the  DOA  estimation 
process.  The  algorithm  is  derived  as  a  batch  method.  A  sequential  implementation  obtained 
by  stepping  through  the  data  in  short  batches  is  applied  to  acoustic  array  time  series  data 
from  held  tests  conducted  by  the  Army  Research  Laboratory. 


1.  INTRODUCTION 

Acoustic  sensor  arrays  offer  the  potential  to  detect  and  track  enemy  ground  targets  such  as  tanks, 
trucks,  and  other  military  vehicles  [2]-[7].  The  acoustic  emissions  of  ground  vehicles  arc  generally  broad¬ 
band  in  nature  and  arc  composed  of  several  discrete  harmonically  related  components  and  an  underlying 
lower  level  broadband  component.  A  typical  plot  of  the  a  target  frequency  spectrum  vs.  time  as  the 
target  moves  past  the  array  is  shown  in  Figure  1.  Traditional  multitarget  tracking  techniques  partition 


Report  Documentation  Page 


Report  Date  Report  Type 

00  Oct  2001  N/A 

Dates  Covered  (from..,  to) 

Title  and  Subtitle 

Contract  Number 

Maximum  A  Posteriori  (MAP)  Multitarget  Tracking  for 
Broadband  Aeroacoustic  Sources 

Grant  Number 

Program  Element  Number 

Author(s) 

Project  Number 

Task  Number 

Work  Unit  Number 

Performing  Organization  Name(s)  and  Address(es) 

George  Mason  University  Kristine  L.  Bell  Dept,  of  Applied  and 
Engineering  Statistics  Fairfax,  VA  22030-4444 

Performing  Organization  Report  Number 

Sponsoring/Monitoring  Agency  Name(s)  and  Address(es) 

Department  of  the  Army,  CECOM  RDEC  Night  Vision  & 

Sponsor/Monitor’s  Acronym(s) 

Electronic  Sensors  Directorate  AMSEL-RD-NV-D  10221 

Burbeck  Road  Ft.  Belvoir,  VA  22060-5806 

Sponsor/Monitor’s  Report  Number(s) 

Distribution/ Availability  Statement 

Approved  for  public  release,  distribution  unlimited 

Supplementary  Notes 

See  also  ADM201471,  Papers  from  the  Meeting  of  the  MSS  Specialty  Group  on  Battlefield  Acoustic  and  Seismic 
Sensing,  Magnetic  and  Electric  Field  Sensors  (2001)  Held  in  Applied  Physics  Lab,  Johns  Hopkins  Univ,  Laurel,  MD  on 
24-26  Oct  2001.  Volume  2  (Also  includes  1999  and  2000  Meetings),  The  original  document  contains  color  images. 

Abstract 

Subject  Terms 

Report  Classification 

unclassified 

Classification  of  this  page 

unclassified 

Classification  of  Abstract 

unclassified 

Limitation  of  Abstract 

UU 

Number  of  Pages 

15 

Frequency  (Hz) 

Figure  1:  Typical  spectrogram  for  an  aeroacoustic  ground  target. 

the  track  estimation  operation  into  two  isolated  processes:  broadband  direction-of-arrival  (DOA)  estima¬ 
tion  from  array  snapshot  data,  followed  by  track  estimation  from  the  DOA  estimates.  This  partitioning 
results  in  procedures  which  arc  suboptimal  and  require  data  association  to  match  DOA  estimates  to  tar¬ 
gets.  In  this  paper,  a  multiple  target  track  estimation  method  that  operates  directly  from  broadband  array 
data  is  presented.  The  maximum  a-posteriori  (MAP)  estimator  for  contact  states  is  derived  for  tempo¬ 
rally  uncorrelated  broadband  target  signals  and  uncorrelated  target  tracks,  where  the  number  of  targets 
is  assumed  known  and  fixed.  This  technique  extends  the  MAP  multitarget  tracking  approach  developed 
for  narrowband  signals  in  [1],  The  track  estimator  is  an  iterative  algorithm  employing  a  nonlinear  pro¬ 
gramming  (NLP)  penalty  method  in  conjunction  with  an  alternating  maximization  algorithm  for  obtaining 
penalized  maximum  likelihood  (PML)  DOA  estimates.  The  penalty  function  couples  the  DOA  estimates 
from  the  PML  algorithm  to  the  tracker  as  synthetic  measurements,  eliminating  the  data  association  step 
of  traditional  multitarget  tracking  approaches.  It  also  creates  a  feedback  mechanism  to  enhance  the  DOA 
estimation  process. 

This  paper  is  organized  as  follows.  In  Section  2,  the  signal  and  motion  models  are  developed  for 
the  broadband  multitarget  tracking  problem.  The  broadband  MAP  track  estimation  method  is  presented  in 
Section  3.  In  Section  4,  a  sequential  implementation  obtained  by  stepping  through  the  data  in  short  batches 
is  described  and  tested  on  acoustic  array  time  series  data  from  field  tests  conducted  at  Aberdeen  Proving 
Ground  by  the  Acoustic  Signal  Processing  Branch  of  the  Army  Research  Laboratory  (ARL).  Section  5 
contains  a  summary  and  areas  for  future  research. 


2.  STATISTICAL  MODEL  AND  ASSUMPTIONS 


We  consider  the  multitarget  tracking  problem  where  there  arc  M  contacts  radiating  broadband  signals 
received  by  an  array  of  N  sensors.  The  number  of  objects  M  is  assumed  known  and  the  trajectories  of  the 
objects  arc  assumed  to  be  uncorrelated  with  the  trajectories  of  other  objects.  The  targets  and  the  array  arc 
assumed  to  lie  in  the  x  —  y  plane.  The  two-dimensional  target  state  is  defined  as  its  bearing  u  =  cos(9) 
and  bearing  rate  ii.  Thus  the  state  of  the  mth  contact  at  snapshot  k  is  x/^m  =  [u/,vm ,  Uk,m]T-  We  assume 
the  motion  of  the  objects  is  described  by  a  first  order  Gauss-Markov  process,  i.e.  for  the  mth  contact. 
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where  F  = 


1  At 
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At  is  the  time  interval  from  k  to  k  +  1,  and  w/,vm  is  a  zero  mean  white  Gaussian 

noise  process  with  covariance  matrix  Q  which  is  assumed  known  and  fixed  over  the  observation  period 
and  equal  for  all  objects.  Under  these  assumptions,  the  probability  density  function  (pdf)  of  m  given 

xfc,m— 1  IS 


P(xfc,m|xfe,m— l)  — 


®XP  {  2  Exfc,»n— l)  Q  (xfc,m  Exfc,m—  l)} 


2vr|Q 


(2) 


where  |Q|  denotes  the  determinant  of  Q.  There  are  1C  snapshots  in  an  observation  batch.  No  data  is 
available  at  k  =  0  so  we  assume  the  prior  distribution  on  initial  object  states  is  Gaussian  with  mean  x o,m 
and  covariance 
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At  the  array,  we  assume  that  the  data  has  been  transformed  into  the  frequency  domain  and  that  there 
arc  L  frequency  bins  of  interest.  The  N  x  1  observed  data  vector  in  the  Zth  frequency  bin  during  the  A:th 
observation  snapshot  has  the  form 


M 

Yk,l  =  X  Sk,l,m^l(Uk,m)  +  n k,h  (4) 

m= 1 


where  $k.i.m  is  a  random  signal  from  the  mth  source  in  the  Zth  frequency  bin  at  the  A:th  snapshot  with 
E\skj.m$l  i  m]  =  (%k,l,m-  The  vector  v/('u/;.-)rt)  is  the  N  x  1  array  response  vector  for  the  Zth  frequency 
bin  to  the  DOA  Uk}m-  and  n/.-./  is  a  N  x  1  vector  of  uncorrelated  sensor  noise  samples  in  the  Zth  frequency 
bin  at  the  Ath  snapshot.  The  source  signals  and  noise  are  assumed  to  be  sample  functions  of  independent 
zero-mean  Gaussian  random  processes.  It  is  assumed  that  the  observations  arc  independent  from  snapshot 
to  snapshot  and  frequency  bin  to  frequency  bin.  The  signal  powers,  m,  arc  assumed  to  be  unknown 
and  to  vary  across  the  frequency  band  and  in  time.  The  noise  covariance  matrix  is  assumed  to  be  known 
and  vary  across  frequency  with  A’[n/,:,/n////]  =  erf  I.  Note  that  the  array  data  depends  on  the  target  state 
xfc,m  only  through  the  bearing  u/iVm  and  not  the  bearing  rate  We  will  use  the  notation 


Hx 
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with  H  =  [1  0],  to  denote  bearing  component  of  the  state  in  the  subsequent  derivation. 


We  denote  the  collection  of  target  states  across  all  targets  at  time  k  as  x/,.  =  {x/,.| .  x/,.2, . . . ,  x/..j\/} 
and  the  collection  of  target  powers  at  time  k  in  frequency  bin  l  as  ak,i  =  {akii,  a^  i^,  ■  ■  ■ ,  afc,z,M}-  The 
collection  of  tai'get  powers  over  tai'gets  and  frequency  at  time  k  is  defined  as  ak  =  {akj  1,  Oik, 2,  ■  ■  ■ ,  c xk,L }, 
and  the  collection  of  data  vectors  across  frequency  at  time  k  is  defined  as  y/,:  =  {y/,.  | .  y*.^,  . . . ,  y 
It  is  also  useful  to  define  the  collection  of  target  states  over  the  batch  (the  track)  for  each  tai'get  as 
=  {x0,m,Xiim,  .  .  . 


At  each  snapshot  and  in  each  frequency  bin,  the  array  data  ykj  is  then  jointly  complex  Gaussian  with 
zero  mean  and  covariance  matrix 


M 
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and  the  pdf  of  the  array  data  conditioned  on  the  tai'get  states  is  given  by 
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We  have  used  the  notation  p(yk,i  |xfc  :  ak,i)  to  emphasize  that  this  pdf  is  conditioned  on  the  random  vector 
Xfc  but  is  also  a  function  of  the  non-random  but  unknown  parameter  vector  ak,l-  At  each  snapshot  the  joint 
pdf  of  the  broadband  data  conditioned  on  tai'get  states  is  the  product  of  the  pdfs  in  each  frequency  bin  and 
is  given  by: 

L 

p{yk\*k  ■  ak)  =  fjp(yfc,z|xfc  :  ak,i).  (8) 

l=i 

The  single  snapshot  joint  pdf  of  the  observations  and  contact  state  conditioned  on  the  previous  contact 
state  is  then 

M 

p{yk,^k\^k-i  ■  ock)  =  p{yk\*k  ■  OLk)  p(xfc  ,m|xfe;m— l)-  (9) 

m—  1 

Let  X,  A,  and  Y  denote  the  batch  collections  X  =  {Xi,  X2, . . . ,  Xm}  =  {xi,  X2, . . . ,  x^;},  A  = 
{ai,  a-2, ... ,  0.) c},  and  Y  =  {yi, . . . ,  y/c}-  The  joint  pdf  of  the  array  snapshot  data  and  tai'get  states  over 
the  batch  is  given  by 
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3.  BROADBAND  BATCH  MAP  ALGORITHM 

In  the  classical  single  tai'get  tracking  problem  where  the  observations  are  a  linear  function  of  the  tai'get 
states  and  the  observations  and  states  are  Gaussion,  the  discrete  time  Kalman  filter  provides  the  mini¬ 
mum  mean  square  error  (MMSE)  and  MAP  estimates  of  the  tai'get  states  given  the  observations.  When  a 
batch  of  observations  is  used,  the  MMSE  and  MAP  estimates  are  obtained  from  the  fixed  interval  Kalman 
smoother  [8].  In  the  problem  considered  here,  the  observations  depend  on  the  target  states  in  a  nonlinear 


manner,  therefore  the  MMSE  and  MAP  estimates  will  yield  different  solutions.  We  also  have  the  added 
complication  of  the  unknown  nuisance  parameter  vector  A.  The  MAP  methodology  provides  a  tractable 
framework  for  solving  this  problem.  We  can  jointly  find  the  MAP  estimate  of  X  and  the  maximum  like¬ 
lihood  (ML)  estimate  of  A  by  maximizing  the  joint  pdf  p( Y,  X  :  A),  or  equivalently  its  logarithm,  with 
respect  to  both  X  and  A.  In  [1]  a  practical  algorithm  for  batch  track  estimation  for  the  narrowband  prob¬ 
lem  was  developed  using  the  penalty  method  of  nonlinear-  programming.  Earlier  versions  also  appeared  in 
[9]  and  [10].  In  this  section,  a  batch  technique  for  broadband  data  is  presented. 

The  MAP/ML  estimates  of  X  and  A  are  the  solutions  to  the  optimization  problem: 

/  K  M 

•  ( np(yfcixfe :  ^  n 

\k=l  m=  1 

Following  [1],  to  assist  in  the  solution  we  introduce  a  set  of  auxiliary  DOA  variables  fjik,m  for  each  tar¬ 
get  and  snapshot,  and  define  collections  of  these  variables  as  p,k  =  {pk,i,  Pk, 2,  •  •  • ,  Pk,M }  and  M  = 
{/^i,  /r2, . . . ,  AW}-  We  replace  uk^m  =  Hx^m  with  the  new  auxiliary  variable  pk,m  in  the  p(yfc|xfc  :  ak) 
term.  In  order  to  retain  the  original  optimization  problem,  we  then  require  the  new  variable  to  be  equal  to 
the  old  variable,  i.e.  Hk,m  =  m  =  Hx^m.  The  unconstrained  optimization  problem  in  Eq.  (1 1)  can  be 
written  as  an  equivalent  constrained  optimization  problem  as  follows: 

'  M  /  K.  M 
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This  formulation  allows  us  to  use  the  penalty  method  of  NLP  (e.g.  [11],  [12])  for  constrained  optimization 
problems.  It  is  an  iterative  procedure  which  involves  solving  a  sequence  of  easier  unconstrained  opti¬ 
mization  problems.  The  easier  problems  are  related  to  the  original  constrained  problem  by  a  continuous, 
differentiable  penalty  function  which  is  equal  to  zero  in  the  feasible  region  where  the  constraints  are  sat¬ 
isfied,  and  which  is  negative  in  the  infeasible  region.  The  penalty  function  relaxes  the  equality  constraint 
resulting  in  a  problem  which  is  an  approximation  to  the  original  problem.  With  each  iteration,  a  stronger 
penalty  is  imposed  for  infeasibility,  and  the  solution  to  the  unconstrained  problem  converges  to  the  solu¬ 
tion  to  the  original  constrained  problem.  An  overview  of  the  method  and  the  convergence  properties  is 
provided  in  [1].  As  in  [1],  we  choose  the  quadratic  penalty  function. 
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where  erf.  is  a  parameter  that  affects  the  strength  of  the  penalty. 


(14) 


To  enforce  a  more  costly  penalty  at  each  iteration,  a  term  ( cq )  1  scales  the  penalty  function,  where 
q  is  the  iteration  index  and  cqiq  =  1,  2, ...  is  a  positive,  decreasing  sequence  converging  to  zero.  The 


penalized  unconstrained  maximization  problem  is  given  by 
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Expanding  and  rearranging  the  terms  in  Eq.  (15),  we  have 
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Note  that  the  first  term  is  only  a  function  of  M  and  A,  the  second  term  is  only  a  function  of  X,  and  the 
third  term  provides  the  coupling  between  the  parameter  sets. 


First  consider  that  for  a  fixed  M  and  A,  we  can  find  X  by  maximizing  over  the  second  and  third  terms 
in  Eq.  (16).  These  terms  decouple  with  respect  to  the  targets,  therefore  the  problem  reduces  to  solving  M 
separate  track  estimation  problems.  Expanding  the  pdfs,  the  problem  becomes 
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This  problem  has  the  form  of  the  classical  single  source  tracking  problem  with  / ik,m  acting  as  the  noisy 
measurements,  and  cqaf.  m  the  measurement  variance.  The  solution  is  the  fixed  interval  Kalman  smoother 
[8].  The  implementation  consists  of  the  standard  forward  Kalman  filter,  followed  by  a  backward  smooth¬ 
ing  filter.  The  forw  ard  Kalman  filter  is  initialized  with  x0|0  =  xo  and  P0|o  =  ^o-  The  state  estimates  and 
their  error  covariance  matrices  are  computed  sequentially  for  k  =  1 , ,IC  using: 
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Once  the  forward  filtering  pass  has  completed,  the  backward  smoothing  pass  updates  the  state  estimates 
and  covariance  matrices  for  k  =  1C  —  1, . . . ,  0  using: 
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(23) 
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The  final  track  estimate  for  the  mth  target  is  Xfc,m  =  ~x-k\1C->  k  =  1 ,IC. 


Next  consider  that  for  a  fixed  X,  we  can  solve  for  both  M  and  A  by  maximizing  over  the  first  and 
third  terms  in  Eq.  (16).  These  terms  decouple  with  respect  to  the  snapshots,  therefore  this  problem  reduces 
to  solving  /C  separate  multiple  source  broadband  DOA  estimation  problems  of  the  form 
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This  is  a  maximum  penalized  likelihood  (MPL)  estimation  problem.  In  [1],  it  was  solved  iteratively  us¬ 
ing  the  expectation-maximization  (EM)  algorithm  [13].  Here,  we  choose  to  solve  it  iteratively  using  the 
relaxation,  or  alternating  maximization  (AM),  method  because  the  AM  method  usually  converges  faster 
than  the  EM  algorithm.  The  derivation  is  given  in  Appendix  A. 


We  then  alternate  between  estimating  the  penalized  DOAs  and  power  estimates  using  the  DOA  AM 
algorithm,  and  the  M  track  estimates  via  M  fixed  interval  Kalman  smoothers.  As  presented  above,  there 
arc  three  levels  of  iteration:  the  penalty  method  iteration  in  which  the  penalty  parameter  forces  the  solu¬ 
tion  into  the  feasible  region,  the  AM  iteration  between  DOA  estimation  and  track  estimation,  and  the  AM 
iteration  used  in  broadband  DOA  estimation.  Each  of  the  iterations  may  be  performed  until  a  convergence 
criterion  is  satisfied  or  for  a  fixed  number  of  cycles.  The  trade-off  is  algorithm  complexity  versus  esti¬ 
mation  accuracy.  In  this  paper,  we  choose  to  perform  the  two  AM  iterations  only  once  for  each  penalty 
method  iteration,  thus  there  is  only  one  global  iteration  loop. 

The  term  cqo\  m  controls  the  strength  of  the  penalty  in  the  DOA  estimation  stage  and  acts  as  the 
measurement  error  variance  in  the  tracking  stage.  In  [1],  it  was  suggested  to  set  of  m  proportional  to  the 
Cramer-Rao  bound  and  to  let  cq  decrease  exponentially.  Here  we  take  a  simpler  approach  in  specifying 
of  m  and  set  it  inversely  proportional  to  the  mth  target’s  total  estimated  power  in  the  kth  snapshot,  i.e. 

<Jk,m  =  ft  ]  Olk,l,ni^j  ■  (2V) 

An  explicit  pseudo-code  description  of  the  batch  algorithm  using  this  method  for  setting  of  m  and  the 
simplified  iteration  strategy  is  provided  in  Table  1 . 


4.  BROADBAND  SEQUENTIAL  MAP  ALGORITHM 
AND  EXPERIMENTAL  RESULTS 

The  batch  algorithm  provides  an  elegant  solution  to  multitarget  tracking  without  requiring  data  as¬ 
sociation.  However,  many  systems  cannot  wait  while  a  batch  of  data  is  collected  and  processed,  and  we 
would  like  a  real-time  (sequential)  solution.  As  suggested  in  [1],  the  batch  method  can  be  extended  to  a 
sequential  method  by  moving  through  the  data  with  a  shorter  batch  window  of  length  /Cs  and  a  stride  of 
length  A k.  In  the  full  batch  method  JCS  =  /C  and  A&:  =  0,  while  a  fully  sequential  method  would  use 


Table  1:  Batch  broadband  MAP  multitarget  tracking  algorithm  pseudo  code. 
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Track  Estimates:  Fixed  Interval  Kalman  Smoother 
for  m  =  1 , . . . ,  M 

Initialize  x0|0  =  x0,m,  P0|0  =  ^o,m 
for  k  =  1 , . . . ,  K, 
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end  {m} 
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JCS  =  1  and  A k  =  1.  In  between,  we  have  batch-sequential  methods  where  A k  is  determined  by  how  of¬ 
ten  state  updates  arc  needed,  and  K,s  is  chosen  to  balance  estimation  accuracy  with  algorithmic  complexity. 


In  the  ARL  data  sets,  most  of  the  acoustic  energy  of  the  targets  is  concentrated  in  the  frequency  band 
below  200  Hz.  The  ground  vehicles  exhibit  both  broadband  energy  and  strong  harmonics  which  arc  non¬ 
stationary  due  to  vehicle  maneuvering  and  environmental  factors.  The  sensor  array  is  a  seven  element 
circular  array  consisting  of  six  elements  in  a  circle  of  radius  4  ft.,  plus  one  sensor  at  the  center.  The  array 
is  placed  very  close  to  the  road  traveled  by  the  vehicles,  therefore  the  range  and  bearing  of  the  targets 
changes  rapidly  as  the  target  passes  by  the  sensor.  Furthermore,  the  signal  power  level  varies  with  range, 
thus  the  multiple  target  scenarios  have  sources  with  highly  variable  power  levels  as  well  as  fluctuating 
DOAs.  The  data  is  sampled  at  a  frequency  of  1024  Hz.  The  1024  samples  are  converted  to  the  frequency 
domain  using  a  1024-point  FFT  to  provide  adequate  frequency  resolution,  resulting  in  one  frequency  do¬ 
main  data  snapshot  per  second.  Previous  studies  have  shown  that  reasonably  good  DOA  estimates  can 
be  obtained  using  the  data  from  40-85  Hz.  Using  higher  frequency  data  provides  increased  resolution  at 
the  expense  of  global  estimation  errors  due  to  grating  lobes  in  the  array  beampattern  [4], [5].  The  large 
error  estimates  arc  extremely  detrimental  to  tracking  performance,  thus  we  restrict  processing  to  the  lower 
frequency  bins  to  avoid  these  errors.  The  noise  power  (which  is  assumed  known)  was  estimated  during 
from  a  quiet  interval  in  the  data  when  no  targets  were  present. 

The  data  was  processed  in  blocks  of  Ks  =  60  seconds  with  a  stride  of  A k  =  10  seconds.  The  block 
length  was  chosen  to  provide  sufficient  track  estimates  for  smoothing  over  periods  when  a  weaker  target 
is  masked  by  a  stronger  target.  The  stride  was  chosen  to  trade  off  providing  track  updates  in  a  timely 
manner  vs.  computational  complexity.  The  blocks  overlapped  by  50  snapshots.  State  estimates  for  the 
overlapping  snapshots  obtained  from  the  previous  block  were  used  as  the  initial  values  for  the  current 
block.  To  initialize  the  state  estimates  for  the  new  snapshots,  we  simply  projected  out  in  time  with  the 
motion  model  from  the  most  recent  state  estimate.  The  batch  algorithm  was  applied  to  the  current  block 
with  qmax  =  2  iterations.  The  final  state  estimates  from  this  block  replaced  the  estimates  from  the  previ¬ 
ous  block.  To  ensure  continuity  of  the  final  track  estimates  as  the  block  moved  through  the  data,  the  prior 
mean  and  covariance  matrix  of  the  target  state  for  the  current  block  was  set  equal  to  the  track  estimate  at 
the  snapshot  just  prior  to  the  current  block,  with  a  zero  variance  on  the  DOA  and  a  small  variance  on  the 
DOA  rate.  To  initialize  the  algorithm,  the  hacks  in  the  first  block  were  initialized  to  have  constant  DOA 
and  power  spectrum  equal  to  the  true  value  at  the  beginning  of  the  block. 

Figures  2  and  3  show  results  obtained  for  two  scenarios  involving  two  targets.  Each  figure  shows 
the  true  DOAs  of  the  targets,  the  DOA  estimates  obtained  using  incoherent  minimum  variance  (IMV) 
beamforming  as  in  [4],  the  DOA  estimates  obtained  using  broadband  ML,  and  the  target  tracks  from  the 
broadband  MAP  tracker.  In  the  first  scenario,  the  targets  approach  from  different  directions  and  pass  each 
other  close  to  the  array.  The  targets  arc  widely  separated  in  angle  except  for  a  short  period  of  time  while 
they  arc  crossing,  and  have  nearly  the  same  power  levels.  For  this  scenario,  the  single  snapshot  DOA 
estimation  techniques  perform  quite  well  except  at  the  beginning  of  the  experiment  and  for  a  short  period 
while  the  targets  arc  crossing.  ML  seems  to  give  better  estimates  than  IMV.  The  tracker  is  able  smooth 
out  the  conflicting  data  and  follow  both  contacts  quite  well.  The  tracks  deviate  a  little  from  the  true  DOAs 
during  the  crossing  period,  but  arc  able  to  recover  quickly. 
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Figure  3:  Wideband  DOA  and  MAP  tracker  estimates  for  two  target  following  scenario. 


In  the  second  scenario,  the  targets  follow  one  behind  another.  The  targets  start  out  with  nearly  the  same 
DOA.  As  the  first  target  passes  the  array,  it’s  DOA  changes  rapidly  from  about  270°  to  about  90°  and  it 
is  so  loud  that  it  masks  the  second  target.  The  IMV  technique  does  not  report  any  DOA  estimates  for  the 
weaker  target  during  this  time  and  the  ML  technique  (which  is  required  to  report  two  estimates)  provides 
random  estimates  initially,  then  estimates  which  are  due  to  a  grating  lobe  from  the  stronger  target.  After 
the  first  target  moves  away,  the  targets  arc  widely  separated  and  both  can  be  heard.  Both  DOA  estimators 
work  well.  As  the  second  target  approaches  the  array,  it  masks  the  first  target.  The  IMV  technique  can  no 
longer  detect  the  first  target.  The  ML  technique  again  provides  erroneous  estimates  which  can  be  seen  as 
belonging  to  energy  from  a  grating  lobe  of  the  loud  target,  then  is  able  to  provide  reasonably  good  esti¬ 
mates  of  the  weaker  source.  This  is  a  difficult  scenario  for  the  tracker.  Initially,  the  tracker  cannot  track  the 
weaker  target,  but  quickly  locks  on  after  the  first  target  has  passed.  The  tracker,  like  the  DOA  estimators, 
is  able  to  provide  an  accurate  track  of  the  second  target  as  it  moves  past  the  array.  The  tracker  becomes 
confused  when  the  first  target  is  masked,  but  is  eventually  able  to  recover  and  track  the  first  contact  again. 
Although  the  tracker  is  clearly  influenced  by  the  poor  ML  estimates,  the  penalty  function  in  the  tracker 
DOA  estimator  and  the  smoothing  operation  allow  the  tracker  to  recover. 

These  results  demonstrate  that  the  extra  processing  performed  by  the  tracker  can  improve  detection  and 
localization  of  multiple  sources  over  wideband  DOA  estimation  alone.  The  improvement  in  performance 
is  due  to  both  the  extra  smoothing  provided  by  the  Kalman  smoother,  and  the  focused  DOA  estimation 
provided  by  the  MPL  formulation. 


5.  SUMMARY 

The  narrowband  MAP  multitarget  tracker  developed  in  [1]  was  extended  to  handle  broadband  targets. 
To  construct  this  algorithm,  we  introduced  an  auxiliary  DOA  parameter  and  a  NLP  penalty  technique  to 
decouple  the  Gauss-Markov  motion  model  and  the  array  data  model.  Convergence  of  the  artificial  prob¬ 
lem  to  the  original  problem  was  enforced  through  a  stronger  penalty  with  each  iteration.  In  decoupling 
the  problem,  DOA  and  power  estimates  were  obtained  from  an  AM  algorithm  solution  to  a  maximum  pe¬ 
nalized  likelihood  problem,  and  track  estimates  were  obtained  from  M  single  target  fixed  interval  Kalman 
smoothers,  with  the  DOA  estimates  serving  as  noisy  measurements.  The  batch  algorithm  was  imple¬ 
mented  in  a  sequential  manner  by  stepping  through  the  data  with  short  overlapping  batches. 

Results  from  acoustic  time  series  data  collected  by  ARL  show  that  the  tracker  can  improve  detection 
and  localization  of  multiple  sources  over  wideband  DOA  estimation  alone.  However,  the  tracker  still  has 
trouble  during  periods  when  one  target  is  so  loud  that  it  masks  the  other  target.  One  possible  solution  is 
to  place  the  sensor  array  so  that  it  is  not  so  close  to  roadways  or  routes  taken  by  the  vehicles.  This  would 
prevent  the  array  from  being  overwhelmed  by  one  source,  and  would  decrease  the  rapidity  with  which  the 
DOA  varies  as  it  moves  past  the  array.  Another  possibility  is  to  jointly  process  data  from  multiple  sensor 
arrays  [6], [7],  Not  only  would  this  increase  the  amount  of  data  available,  but  during  periods  when  a  target 
is  close  to  one  array,  it  would  not  be  close  to  the  rest  of  the  arrays,  so  target  masking  will  not  occur  at  all 
arrays  simultaneously. 


APPENDIX  A.  BROADBAND  ALTERNATING  MAXIMIZATION  ALGORITHM 

FOR  MPL  DOA  ESTIMATION 


We  start  from  the  single  snapshot  penalized  likelihood  function  in  Eq.  (26).  Expanding  terms,  it  has 
the  form 


A  =  lnp(yk  :  /rfc,  ak)  =  ^ln 
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To  reduce  notation,  we  let  Kyfc  ;  =  Kyfe  4  (/**,  ctk,l)  and  write  (28)  as 
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We  can  use  the  alternating  maximization  technique  to  iteratively  solve  for  parameters  associated  with 
one  of  the  M  targets  by  holding  all  of  the  parameters  associated  with  the  other  targets  fixed.  The  likelihood 
increases  at  each  iteration  and  convergence  to  a  local  maximum  is  guaranteed.  This  reduces  to  a  series 
of  MPL  problems  involving  a  single  source  in  known  colored  Gaussian  noise.  Consider  estimation  of  the 


mth  target’s  DOA  Hk.rn 

and  power  spectrum  ce^,i ••• ,  &k,L,m-  We  ca  rewrite  Kyfc  t  as 
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m'^m 

is  the  colored  noise  covariance  matrix.  It  depends  only  on  the  noise  power  and  the  parameters  of  the  other 
sources,  which  arc  assumed  known.  Performing  a  pre-whitening  operation  on  the  data,  we  have 


y  k,l  =  E/J/riW  (32) 


The  pre-whitened  data  covariance  matrix  becomes 
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where 

The  penalized  log-likelihood  function  in  (29)  can  be  written  as 

A  =  E  {-y&K^yfc,*}  -  Eln  \K9kA  -  {fik,m2c  f* 
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where  7  is  a  term  which  does  not  depend  on  the  parameters  of  the  mth  target.  Using  the  definition  in  (33) 
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and  the  penalized  log-likelihood  reduces  to  (dropping  unnecessary  terms): 
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Maximizing  over  ak,i,m  first,  we  take  the  derivative  and  obtain 
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To  reduce  spurious  estimates  due  to  false  alarms,  we  restrict  the  value  of  the  signal  power  estimate  to  be 
no  less  than  the  noise  power.  Our  constrained  estimate  becomes: 

1 2 


6(k,l,m{Hk,m)  —  max  (7 1 


1 


( 


I  v;(Wfc,m)r 


y  k,M»k,r. 


V 


I V;  (/Tfc,m)  | 


-  1 


(41) 


To  find  an  estimate  of  we  substitute  the  estimates  ak,i,m{Hk,m)  for  l  =  1  back  into  the 

penalized  log  likelihood  function  in  (38)  and  maximize  with  respect  to  / ik,m •  There  is  no  closed  form 
expression,  so  a  one-dimensional  search  is  required. 


f^k,m 


argmax 


E  |y  fcW/9c,m)r 


&k,l,m(lJjk,r 


1=1 


^ln{l  +&k,i  ,m  (M k  ,ra  )  I  J* 


1  +  &k,l,m{l^k,m)  |  V/(/i/c,m)  |  J 
2 1  (/^/c,m  H  X/^m) 


1=1 


^ CQak,m 


(42) 


The  power  estimates  arc  then  found  by  substituting  the  DOA  estimate  fik,m  back  into  (41)  for  /  =  1 , . . . ,  L. 
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